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Abstract — Given two independent point processes and a certain 
rule for matctiing points between them, what is the fraction of 
matched points over infinitely long streams? In many application 
contexts, e.g., secure networking, a meaningful matching rule is 
that of a maximum causal delay, and the problem is related to 
embedding a flow of packets in cover traffic such that no traffic 
analysis can detect it. We study the best undetectable embedding 
policy and the corresponding maximum flow rate — that we call 
the embedding capacity — under the assumption that the cover 
traffic can be modeled as arbitrary renewal processes. We find 
that computing the embedding capacity requires the inversion 
of very structured linear systems that, for a broad range of 
renewal models encountered in practice, admits a fully analytical 
expression in terms of the renewal function of the processes. 
Our main theoretical contribution is a simple closed form of 
such relationship. This result enables us to explore properties 
of the embedding capacity, obtaining closed-form solutions for 
selected distribution families and a suite of sufficient conditions 
on the capacity ordering. We evaluate our solution on real 
network traces, which shows a noticeable match for tight 
delay constraints. A gap between the predicted and the actual 
embedding capacities appears for looser constraints, and further 
investigation reveals that it is caused by inaccuracy of the 
renewal traffic model rather than of the solution itself. 



I. Introduction 

CONSIDER the pair of timing sequences represented by 
the point processes S and T in Fig. [H where points 
are matched according to some prescribed rule. What is the 
maximum achievable fraction of matched points (embedding 
capacity) given the two processes and the matching rule? 
How do statistical properties of the point processes affect the 
maximum fraction of matching? The main theme of this paper 
is that of providing analytical tools for computing the embed- 
ding capacity of two independent and identically distributed 
renewal processes, when the coupling rule is formulated in 
terms of a causal delay constraint. 

The above problem naturally arises in many applicative 
scenarios: from intelligence applications aimed at tracing 
relationships among individuals (e.g., in social networks), to 
the discovering of neuron connections by measurements of 
firing sequences, and so forth lH], An application closer 
to the communication area concerns the anonymous relaying 
of messages in distributed architectures, or the detection of 
clandestine information flows in wireless systems. In fact, 
the evaluation of the embedding capacity under causal delay 
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Fig. 1. Notional sketch of the addressed problem, with arrival epochs of 
processes S and T matched according to a delay constraint A. Matched 
points are marked by circles, unmatched by diamonds. 



constraint has been recognized as a relevant problem in the 
context of secure networking, where the focus is on informa- 
tion flowing that is anonymous with respect to an attacking 
eavesdropper [3J, or, in a reversed perspective, clandestine 
with respect to a legitimate traffic analyst 13. 

In these contexts — to which we specifically refer in the 
paper — the two processes represent the sequences of time 
epochs (traffic patterns) at which successive packets leave two 
nodes of the network and, for security requirements, packets 
are encrypted so that they do not reveal special characteristics. 
Still, the act of transmission itself cannot be kept secret, and 
timing analysis can be performed. 

Given that nodes are unable to hide the act of transmission, 
they must hide the information flow packets into their normal 
transmission scheduling, which provide cover traffic for the 
desired flow. The nodes can mask the timing relationships 
by properly delaying the transmission of information packets 
and/or multiplexing information packets with dummy packets 
or packets from other flows. With a sufficient amount of 
perturbation, an information flow can be disguised as traffic of 
arbitrary patterns. In particular, the flow can appear identical to 
independent traffic following certain transmission schedules. 

As a consequence, every transmission schedule (or cover 
traffic) has certain capacity of being utilized to transmit infor- 
mation flows covertly. The matching capability of a particular 
schedule takes the operational meaning of an embedding 
capacity, that is, the maximum fraction of information packets 
that can be embedded in the cover traffic following this 
schedule, leaving no chances of discovering the presence of 
the flow itself. 
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A. Summary of Results 

The embedding capacity for a Poisson process under causal 
delay constraint is known, see [4|. The Poisson assumption, 
however, rarely fits real traffic and, to date, analytical formulas 
for arbitrary renewal traffic are still missing. The contribution 
of this paper is in filling this gap. 

We find that the embedding capacity is related to the 
invariant distribution of a certain Markov chain. First, we 
prove the existence of such distribution, so that capacity 
evaluation requires the solution of an integral equation. We 
attack this problem by exploiting the powerful tools offered 
by the Riemann-Hilbert theory, which allows us to derive the 
following approximation for the embedding capacity; 
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where A is the rate of the processes, A is the delay constraint, 
and m{t) is the renewal function of the (scaled to unit rate) 
underlying process. The accuracy of this formula is excellent 
for a very broad range of renewal processes of interest for 
the applications, see Sect. |V] We also show how C* can be 
computed to any degree of approximation by inverting a very 
structured linear system, and provide a first-order correction 
expressed in closed form. 

The significance of the above formula is that C* depends 
only on the renewal function which is the key quantity in 
renewal theory, as such, is well studied and understood. In 
many cases of practical interest, the integral involved can 
be evaluated explicitly, from which physical insights can be 
gained. 

The above expression is then used to relate the physical 
parameters and properties of the renewals to the embedding 
performance. In the asymptotic regime of large AA, the dis- 
persion index 7 is the only relevant quantity, and the capacity 
scales as 1 — 7/(AA). Stochastic variability is instead the key 
(for any AA) to compare different interarrival distributions: 
less variable interarrivals yield a larger embedding capacity. 

B. Relevance to Secure Networking 

One applicative scenario of interest is that of secure net- 
working. Consider two packet streams in a network, whose 
transmission timestamps are represented by point processes 
S = (S'i,52,...) and T = {Ti,T2,...). We assume that 
the packet content is fully protected by encryption, while 
the transmission patterns are relatively easy to obtain by 
a monitoring agent, which is capable of performing traffic 
(actually timing) analysis. 

In one possible scenario, S and T are transmission activities 
of two nodes A'^i and N2 in a wireless network. The existence 
of a flow from 5 to T implies that N2 is acting as a relay for 
(part of the transmissions from) A^i, thus revealing a multi- 
hop route that is otherwise unobservable in protocol or content 
domain. The monitoring agent is interested in discovering 
whether or not such relaying exists; conversely, the nodes 
would like to hide the presence of the information flow (e.g., 
to preserve anonymity) by transmitting independently. 



In another case, S represents the pattern of the packets 
transmitted to a multiaccess relay through an ingoing link Li, 
and T is the pattern for an outgoing link L2, both among 
multiple ingoing/outgoing links. It is known that information 
is flowing from Li through the multiaccess relay, but it is not 
known whether or not the outgoing link L2 is being used for 
this. The monitoring agent is interested in tracing the route of 
the flow, and the multiaccess relay intends to hide the route 
by scheduling the two links independently. 

C. Related Work & Organization 

The roots of packet embedding into cover traffic can be 
traced back to the early '80s. The problem of avoiding traffic 
analysis using special relay policies was first considered in |5 1, 
with the adoption of the so-called MIX relays, that perform 
multiplexing, scrambling and encryption of the incoming 
traffic in order to eliminate the correlation with the outgoing 
traffic. Since then, several studies have been made in order to 
improve relay performances, see e.g. ||6], (|7]. More recently, 
it has been shown how statistically independent transmission 
schedules can achieve perfectly anonymous relaying, with 
emphasis on the maximization of the carried information 
capacity (3). 

Also related to our problem is the network security issue 
referred to as stepping-stone attack B, ||9], in which an adver- 
sary launches an attack through a sequence of compromised 
servers, and one would like to trace the sequence to the origin 
of the attack. For wireless networks, an ad hoc network may 
be subject to the worm-hole attack [TOl, where the attacker 
hijacks the packets of a node and channels them through a 
covert tunnel. In such scenarios, the maximum information 
rate sustainable by the attackers is related to the embedding 
capacity of the node traffic patterns. 

From an information theoretic perspective, the problem of 
secure communications, in terms of maximizing the reliable 
rate to a legitimate receiver with secrecy constraints with 
respect to an eavesdropper, has been extensively studied, since 
the pioneering works (TT), IIT2I . |fT3l , up to recent extensions, 
including multiaccess [14j, fading |[T5| , feedback |,16J , and 
broadcast ITtI channels, among many others. We stress that 
the specific scenario of interest for this paper is instead secure 
networking with focus on anonymous relaying of information, 
according to the model proposed in ||3], IS- 

Formal studies of the embedding properties of renewals have 
been carried out in [|3], (|4] , with extensions to distributed de- 
tection with communication constraints IITSl , |[T9l . The authors 
of f4l settled up the problem from the traffic analyzer's per- 
spective, where the role of the embedding capacity is replaced 
by that of undetectable flow. They found a closed formula for 
the capacity under the Poisson regime. General renewal traffic 
models in many applications (inside the communication area 
as well as outside that) are far from being approximated as 
Poisson, such that several extensions of the above studies in 
this direction have been proposed, see ||20) , 11211 . However, a 
tractable analytical formula for the embedding capacity under 
arbitrary renewal traffic is still missing. 

The remainder of this paper is organized as follows. Sec- 
tion [n] formalizes the problem, the main results of the paper 
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are presented in Sect. |III] and Sect. |IV]is devoted to the main 
mathematical derivations. Sect. |V] concerns the appUcation 
of the main theoretical findings to specific examples, while 
Sect. |VT] addresses the problem of classification and ordering 
of renewal processes in terms of their embedding capacity. 
Finally, Sect. IVIII presents the results of experiments on real 
network traces, and conclusions follow in Sect. IVIIII An 
appendix contains some mathematical derivations. 

II. Problem Formulation 

Capital letters denote random variables, and the correspond- 
ing lowercase the associate realizations, while Pr and E denote 
probability and expectation operators, respectively. 

Consider two point processes S — {Si, S2, • • ■ ) and T ~ 
(Ti, T2, . . . ) defined over the semi-axis t e (0, 00). Points that 
are matched over the two processes form an information flow 
in the sense that one point in a matched pair can be thought 
of as a relayed copy of the other. We are interested in delay- 
sensitive directional flows, for which matched points obey a 
causal bounded delay constraint as follows. 

Definition 1 (Information flow) Processes 
W = (Wi,W2,...) and Z = (^1,^2,...) fonn a A- 
bounded-delay information flow in the direction W — > Z 
if for every realization {wi\ and {zi\, there is a one-one 
mapping {wi\ — {zj} that satisfies the causal bounded delay 
constraint < — < A, Vi. o 
Here A > is a known constant representing the maximum 
tolerable delay during relaying. 

Given point processes S — (^i, 5*2, . . . ) and T = 
(Ti, T2, . . . ), an information flow can be generated by finding, 
for each realization of the processes, subsequences that admit 
a valid one-one mapping. This is controlled by an embedding 
policy. 

Definition 2 (Embedding policy) An embedding policy e 
selects subsequences W"^ of S and Z'^ of T to form an 
information flow. o 
The name "embedding" is due to the fact that to an outsider 
who cannot observe the selection, it is not known which points 
belong to an information flow or even if there is a flow, and 
thus the flow is embedded in the overall processes (S, T). For 
the same reason, (5, T) is called cover traffic. 

Let £ = {e} be the set of admissible embedding policies. 
Given e G £, the cover traffic {S, T) is decomposed into 

S = W'®U\ T--Z'®V\ 

where (W^, Z^) forms a valid information flow. Here ® is the 
superposition operator for point processes: {ci}= {ui} ® {bi} 
means that {ci} — {a^} U {hi} with ci < C2 < . . . . 

Given the cover traffic, each embedding policy has a certain 
capability of hosting information flows, quantified as follows. 

Definition 3 (Efficiency) Given cover traffic (5, T), the 
efficiency of an embedding policy e G £ is measured by 



where A^w" (0> ^Z' (^) '^'"^ ^'^^ counting processes for the 
embedded information flow, so are Ng{t), N-j-{t) for the cover 
traffic (assuming the limit exists almost surely). o 
That is, the efficiency is the asymptotic fraction of matched 
points in the cover traffic. We are interested in the highest 
efficiency that we call the embedding capacity. 



Definition 4 (Embedding capacity) C* 
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The embedding capacity C* is a function of the cover traffic 
and the flow constraints (e.g.. A), omitted for simplicity. 
We shall focus on the case that the cover traffic S and 
T are independent and identically distributed (i.i.d.) renewal 
processes, with interarrival random variables X and Y, re- 
spectively. Throughout the paper it is assumed that X and 
Y are absolutely continuous with known Probability Density 
Function (PDF) f{t) and Cumulative Distribution Function 
(CDF) F{t), and that the rate of the processes, denoted by A, 
is finite and nonzero, i.e., < A = 1/E[X] 1/E[y] < 00. 
When the second moment is finite, we define the dispersion 
index as 



7 = A^ VAR[X] = A^ VAR[F] < 00. 



(1) 



ri(e) 



t"-^S3 Ns{t)+Nr{t) 



lim 



III. Characterization of the Embedding capacity 

A. Optimal Embedding Policy 

As a first step toward embedding capacity evaluation, we 
need to find an optimal embedding policy that maximizes the 
number of matched points for any given cover traffic, thus 
achieving the embedding capacity. This has been achieved 
by an existing algorithm called the Bounded Greedy Match 
(BGM) \22\ . It is a simple algorithm that classifies the points of 
two arbitrary point processes as "matched" and "unmatched" 
by sequentially matching points in the two processes under 
a causal delay constraint A and marking the points violating 
this constraint as unmatched. 

The BGM algorithm works as follows. Given realizations 
of two point processes, all the points initially "undetermined", 
the BGM repeats the following steps, see Fig. [T] 

1) Consider the first (in the direction of increasing time) 
undetermined point in the first process, say p^^^; 

2) Find the first undetermined point in the second process 
in the interval [p*-^-*, p^^-* + A], if any, denoted by p'^-*; 

3) If such a point exists, mark both p^^^ and p'^^ as 
"matched"; otherwise, mark p'^^ as "unmatched"; in 
either case, mark all undetermined points in the second 
process before p^^^^ "unmatched". 

Matched and unmatched points are also referred to as "flow" 
and "chaff", respectively. The BGM is optimal in the sense 
that, given two arbitrary realizations of point processes and an 
arbitrary value of A, the algorithm finds the maximum number 
of matched points satisfying the delay bound lE], BJ, ll22l . 

B. Embedding capacity in terms of a Markov Chain 

Our second step in deriving the embedding capacity C* 
consists of modeling the behavior of BGM by a Markov chain, 
whose stationary distribution is directly related to C* . 
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Fig. 2. Three situations arising from applying the BGM procedure to point processes S and T. Chaff points are denoted by "o". Left: the point at Si is 
unmatched, and it is a chaff point in the first process. Center: aU points are matched (no chaff). Right: a chaff point is present in the second process. 



With reference to Fig. |2] let us consider the time difference 
between the first points, that is Zi = Ti — Si. According to 
the BGM algorithm, we have the following three possibilities: 
(i) If Zi > A, the points cannot be matched, and the one in 
S is labeled as chaff. To decide the nature (chaff/non chaff) 
of the point in T, we must check whether it can be matched 
to the next arrival in S, thus computing, see Fig. |2l^a), 

Z2 = Ti — S2 = Zi — X, 

where X is the random variable representing the interarrivals 
in S. 

(a) If < Zi < A, the points match. To check the nature 
of the next incoming points, we update the process as, see 
Fig. life), 

Z2=T2-S2 = Zi+Y-X, 

where Y is the random variable representing the interarrivals 
in T (recall that Y has the same distribution as X). 
(in) If Zi < 0, the points cannot be matched, and the one in 
T is labeled as chaff. To decide the nature of point in S, we 
must check whether it can be matched to the next arrival in 
T, thus computing, see Fig. |2jc). 



Z2 — T2 — Si — Z\ 



Y. 



By repeating for the successive points, we see that a Markov 
process can be compactly defined in terms of the original 
renewals by the following recursion rule 



Zn-1 
Zn-1 



Yn , 



if Z„_i > A, 

if < Zn-1 < A, (2) 

if < 0, 



where Xn and Yn are the interarrivals of the first and the 
second process at the nth step of the chain, following the 
common PDF f{t). 

The Markov chain defined in ^ is schematically illustrated 
in Fig. [3] According to the constitutive equation dU, the 



increment of the Markov chain is the interarrival difference 
Y—X if the chain is currently between and A, which implies 
the matching is successful (e.g., Zi, Zg); the increment is Y 
if the chain is below 0, in which case the reference point 
in the second process is marked as chaff and the reference 
point in the first process remains the same (e.g., Z3, Z4); 
similarly, the increment is —X if the chain is above A, when 
the reference point in the first process becomes chaff and that 
in the second process remains unchanged (e.g., Z7). Note that 
the number of steps of the Markov chain lying inside (resp. 
outside) the barriers and A defines the number of flow (resp. 
chaff) points marked by the BGM algorithm. This suggests 
that a relationship exists between the asymptotic distribution 
of the chain and the fraction of flow points, i.e., the embedding 
capacity. This is made precise in the next section. 



C. Main Results 

The first theorem we present, whose proof is deferred to 
appendix |A] establishes a connection between the embedding 
capacity and the invariant distribution of the BGM Markov 
chain, expressed as the solution of an integral equation. 

Theorem 1 (C* by Markov chain) Let S and T be two 
independent and identically distributed renewal processes, 
with interarrival PDF f{t). Let A be the delay constraint, 
and define a Markov chain by (O. Assume BGM can match 
at least one pair of points in S and T almost surely. 

a) The invariant PDF h{t) of the Markov chain exists and 
solves the following homogeneous Fredholm integral equation 
of the second kind i l25l/ 

/O r+00 
h{T)f{t - T)dT + / h{T)fiT - t)dT 

-00 J A 



h{T)fa{t-T)di 



(3) 
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where fo{t) is the convolution between f{t) and f{-~t). 



defined as f{T).f{T ~ t)dT. 

h) The embedding capacity can be written as 



C* 



l + J^^h{t)dt 



(4) 





Since now, we shall assume that the hypotheses of Theo- 
rem 1 are in force. The next theorems, whose proofs are given 
in Sect. IIVI accordingly focus on the solution of the integral 
equation relevant to capacity computation. We first introduce 
the following definitions. 

Definition 5 (u-PDF) The probability density function k{t) 
of the interarrivals scaled to unit mean, that is, the random 
variables XX and XY, will be called u-PDF. o 

Definition 6 (u-RF) The Renewal Function 

m{t) :=E[iV(i)], 

where N{t) is the number of arrivals in (0, t) of the processes 
scaled to unit rate, having interarrival random variables \X 
and \Y, will be called u-RF. o 

Theorem 2 (Exact value of C*) Under the assumption of 
finite second moment for the interarrivals, the embedding ca- 
pacity of two independent and identically distributed renewal 
processes with rate X, under delay constraint A, is 



C* 



217(0) 



1 + 17(0) 
where r2(/) is the solution of 



(5) 



K{v) 



1 - K{v) 



AAsinc[AA(/ - i/)]di/ 



= AAsinc(AA /) 



1 - 0(0) 



(6) 



K{f) being the Fourier transform of the u-PDF k{t), and 
sinc(t) ~ sm{TT t) / {nt) . 

As a check, let us specialize the above equation to the 
case of exponential interarrivals, for which embedding ca- 
pacity is available in closed form [4]. It is easily seen that 
3? I i^K^f) } ^ allowing direct solution of eq. (|6]l, and 
computation of 11(0) = AA/(2+AA). Substituting into eq. (|5]l, 
this yields 



C* 



AA 



1 + AA 



(exponential), 



that matches the known result from [4 |. 

Note that Theorem 2 still gives an implict solution to the 
problem in terms of an integral equation, which does not 
have a closed-form solution in general. On the other hand, 
eq. (|6]l turns out to be amenable to approximate solutions, 
thus yielding the results stated in the next two theorems. 

Theorem 3 (Approximation of C*} Under the assumption 
of finite second moment for the interarrivals, the embedding 




S 



r 



-Z,^-{Z,+Y) 



Zg- Zj-X 




Fig. 3. Construction of a sample path of the Markov process (lower panel) 
from a realization of the two point processes (upper). In the upper panel, the 
points marked with "o" are those classified as chaff by the BGM algorithm. 



capacity of two independent and identically distributed re- 
newal processes with rate X, under delay constraint A, can 
be approximated as 

" -. (7) 



c* 



1 + 1^ / m{t)dt 



AA 



TO(t) being the u-RF. {> 

Again, let us apply eq. (|7J in the Poisson regime. The u-RF 
of an exponential random variable is m{t) = t, that inserted 
in (|7]i gives 

AA 



C = 



(exponential), 



1 + AA 

implying that, in this particular case, formula (|7) is exact, i.e., 
C* = C. This can be understood by looking at the technique 
used to get the approximatioi{3 in the proof of Theorem 3. 

The relevance of Theorem 3 stems from the fact that, for the 
typical interarrival distributions encountered in many applica- 
tions, the accuracy of the fully analytical approximation (|7]i 
seems to be excellent irrespective of the range of the product 
AA, the tailweight of the distribution, its variance (and even 
for infinite second moment), as confirmed by the examples in 
Sect. [V] Accordingly, Theorem 3 provides us with an accurate 
and mathematically tractable expression for the embedding 
capacity under arbitrary renewal traffic. 

We would like to emphasize that the characterization (|7]l 
relates the sought capacity to the u-RF of the underlying 
process. This highlights the role of the renewal function 
m{t), and reveals that its average ^ m{t)dt is the key 
quantity in determining C. Thus, different traffic models can 
be classified with respect to their embedding capabilities just 
in terms of that average. 

'The terms fc ^ neglected in eq. )23t . are rigorously zero in the 
exponential case, since the integral in l IlOl is zero. 
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We now state a corollary characterizing the asymptotic be- 
havior of the capacity in the limit of A ^ 1/ A. From a known 
property of the renewal function [24], m{t) — t ^ (7 — l)/2 
in the limit of i 00, where 7 is the dispersion index defined 
in eq. ([T]i. Simply plugging that expression in eq. (|7]i would 
give 1 — C ^ 7/(AA). Indeed, we have the following result. 

Corollary 1 (Scaling law for C) Under the assumption 
of finite second moment for the interarrivals, limAA^-oo[l — 
C](AA) = 7, i.e., the embedding capacity in Theorem 3 scales 
as 

7 

AA' 



1-C 





The corollary reveals that, for large values of the product 
AA, the key quantity in determining the capacity is the 
dispersion index: given AA ^ 1, the ability for a type of 
(renewal) traffic to hide information flows in independent 
realizations only depends on the value of the dispersion index 
7, and different traffic models sharing the same dispersion 
index behave similarly. 

Finally, to improve on the approximation in Theorem 3, we 
provide the following theorem that expresses the embedding 
capacity as the solution to a simple linear system. Consider, 
for any integer > 1, the following system 



^ , k\ AA ^ 



k=-N 



AA 



h = -N, ...,N, 



where Ih — ^ for h — 0, and Ih — otherwise. The analytical 
expressions of the entries Ahk, defining a 2N + 1 by 2N + 1 
matrix A, are 



iQO 



i-kk 



Aok = 



A-hk — 



AA 












1 m 
'0 


+ 27rfc 




2(-l)'= 


/ m 
'0 


AA . 



■ AA 



'm{t)dt, 

2TTkt 



AA 

27rfci 
'XA 

(2TTkt\ 



(8) 

dt, k^Q, (9) 
fc^ 0,(10) 



\h — k 



{h-k) 



[M-l)'^ol-*:(-l)'Aot], h^k. (11) 



Theorem 4 (Linear system approximation of C*) Under the 
assumption of finite second moment for the interarrivals, let 
C* — j^^^o) a-^ in Theorem 2. Then, assuming that A 
is invertible, Cl{0) can be approximated as AA/2 times the 
{Q,Q)-entry of matrix A^^, namely $1(0) = ^ {A^^joo- In 
particular, specializing for N = 1, the capacity becomes 



C* 



AA 



2 

AA 



AA 



i{t)dt 



42 

^01 



Am - All 



(12) 







Remark 1. Note that, in the approximation corresponding to 
= 1, a correction term 2 ji^'^°\^^ appears, with respect to 
C in eq. (|7]i, which uses only Aqq. Also, from eqs. dSTl-lfTTTi we 



see that A is very structured and its degrees of freedom grow 
only linearly with N; in fact, A is completely specified by 
assigning one row and the main diagonal. This structure is very 
convenient for numerical tractability. Finally, it is expected 
that the solution becomes more and more accurate as the 
system size N increases. In the section devoted to numerical 
experiments, we show that the zero-order approximation C is 
well satisfying in many cases of interest. Even when this is 
not strictly true, a first-order correction (ITZl offers very good 
results. 

Remark 2. Let us consider a random variable with u-PDF 
k{t) which is zero for t < a, some a > 0. We have m(t) — 
for t < a. This implies that, in the range AA < a, the cross 
terms Aok with k Q vanish, so that the approximation (|7]i 
is exact, and gives the linear relationship C* = AA in the 
considered range, as verified lateiH (see Fig. |6]l. This is also 
consistent with earlier approximation and simulation results 
in mj. 



IV. Proofs of Theorems 2-4 via Riemann-Hilbert 
Theory 

In the following, we make use of suitable normalization 
of the relevant physical quantities, aimed at simplifying the 
mathematical derivation. Indeed, the problem possesses a 
natural scale-free property. For a given distribution of the 
interarrival process, we note that doubling the arrival rate 
"speeds up" the system so that the sample paths can be 
redrawn on a time axis scaled by a factor 2, and halving 
A leaves unchanged the number of matches. We accordingly 
introduce the normalized delay S = AA, and work in terms of 
the unit-mean random variables with u-PDF k{t). It is further 
convenient to symmetrize the problem by shifting the Markov 
chain, as well as the corresponding boundaries, which yields 
to the following Fredholm equation 

u{t) = / u{T)k{t ~ T)dT + / u{T)k{T-t)dT 
J-00 Js/2 
l-S/2 

u{T)ko{t - T)dT, (13) 

-S/2 

where kd{t) is the convolution between k{t) and k{—t). 

Equation ( fT3] l is a homogeneous Fredholm equation of the 
second kind, and we have three different regions where the 
integrals look like convolutions. Were the integral equivalent 
to a convolution as a whole, a simple transform method 
would be directly applicable. To elaborate, let us first con- 
sider what would happen if the function was known within 
the strip [—5/2,5/2]. In this case, the equation would be 
nonhomogeneous, and will be classified as a convolution-type 
equation with two distinct kernels. For this case a powerful 
approach, that can be traced back to Carleman and to Wiener 

-The same conclusion can be also argued as follows. For delay A smaller 
than the minimum allowed interanival time, 5fc — 5fc_i > A, such that the 
probability that matches is the probability that the first aiTival after 5^ 
in T occurs before 5^ -I- A and it can be computed, due to independence 
between the processes, by using the residual lifetime distribution |24|: 
?r[Sk matches] Ri A - F{t)]dt. This implies Pr[Sfe matches] AA, 

for AA < a. By ergodicity, C* = AA in the considered range. 
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and Hopf, still prescribes transforming the equation in the 
Fourier domain [251 • After transformation, the problem falls 
in the class of the so-called Riemann-Hilbert boundary value 
problems. 

The Riemann-Hilbert problem^ in a nutshell, consists in 
finding two functions, analytic in the upper and lower half 
planes, respectively, whose difference on the real axis equals 
a known function (25] f26l. Direct application of this approach 
would require that the sought stationary distribution be known 
within [—6/2,6/2], but this is not our case. Generalizations of 
the method have been proposed. They include the Carleman- 
Vekua regularization method, which suggests to initially treat 
the unknown function as known, and formulating a new 
integral equation in terms of the function in the interval 
[—6/2,6/2], and the work by Jones [27 J, see also [28] . The 
proof that follows is based on these approaches. 

Before, we need some basic notation and concepts about 
one-sided functions and their analytic Fourier transforms, 
which will be useful in the following. For a generic func- 
tion g{t), let g+it) - 9{t)l{t) and g-{t) = -5(t) 
where l{t) is the Heaviside unit-step function l{t) = 1 for 
t > 0, and l(t) = for t < 0. In the Fourier domain, 
this means G+{f ) = J^^°° g{t)e^^''f^dt, and G-(/) = 

By replacing the real parameter / by a complex variable z = 



f+iy, the above integrals become G^{z) = Jq °° g{t)e^'^^^^dt 
and G~{z) = — j'^^g{t)e^'^'^^^dt, which are analytic in those 
regions of the complex plane of the variable z in which they 
are absolutely convergent [25 [: G^{z) is analytic for ^{z) > 
0, and G-{z) for 3(z) < 0. 

From Sokhotski-Plemelj formula [25], or simply decompos- 
ing the Fourier integral into the left and the right part, we have 
that, on the real axis, 

.+ ^,,_Gif) + mGif)] -G(/) + zH[G(/)] 



G+(/) = 



^G-(/) = 



(14) 

where the Hilbert transform 'H[G{f )] = ^ J j^dv has been 
introduced (the integral is in the sense of Cauchy principal 
value). 

A. Proof of Theorem 2 

Consider the the unknown function u{t) in eq. (fTSI l and let 
us define 



where 



L{t) = v+lt - 6/2) -v-{t + 6/2) + uj{t), 



v+{t-6/2) = u{t)l{t- 5/2), 
v-{t + 6/2) = -u{t)l{-t- 6/2), 

u{t) 5/2<t< 5/2 



u{t) 



otherwise 



^Actually, there has been some uncertainty about the original pioneer of the 
approach. According to Musldielishvili L26J "77ie problem formulated above 
is often called the Riemann problem, but the Author considers this name to 
be incorrect [...], because it was first considered by D. Hilbert essentially in 
the form in which it is stated'' 



The corresponding Fourier transforms will be accordingly 
denoted by V+{f), V~{f) and n{f). Note that, from ©, 
we are just interested in h{t)dt — /*^y2 '^{t)dt — J7(0). 

Transforming both sides of the integral equation (fTjt into 
the Fourier domain gives 



V+{f) 



v-{f)e-'^'f + n{f) 



= V+{f)e'''^fK{f) - V-{f)e-'''^fK{f) 

+ ^{m{f)\\ 

where a is the conjugate of a. The above equation can be 
recast as 

l-Kif) ^ l-K{f) 
where we define 



W{f), 



(15) 



1 + 23? 



K{f) 
1 - K{f) 



(16) 



Multiplied by e ^'^^f , eq. (flSl l becomes 



-W{f) 



,-i-iTSf 



l-K{f) l-K{f) 
and using the factorization in (fl4] l: 

yy(j)g-»^<5/ ^ [iy(/)e-»'^*/]+ - [W{f)e'"^f] 
Combining the above equations gives 



V+{f) 



Wif). 



1 - K{!) 

_ V-{f)e-'^'"^^ 
l-K{f) 



[W{f)e 



(17) 



By construction the function ^j^^l-^ is analytic in the upper 
half plane ^{z} > 0, continuous on the real axis, with a 
single pole located at z = 0. By the known property of the 
characteristic function, K'{0) — i2TT, such that this pole has 
order one. On the other hand, the function — i-K{f) — 
analytic in ^{z} < 0, continuous on the real axis, with a single 
pole of order one located at z = 0. Similar considerations 
apply to [W{f)e-'''^^+ and [m/)e-"'^^]-, with the further 
property that there are no poles o 

The asymptotic behavior of the involved functions is essen- 
tially determined by Fourier transforms, such that we assume 
boundedness at infinity. 

Summarizing, the LHS and RHS of eq. ( fTTI ) define functions 
that are analytic in the upper half and lower half planes, 
respectively. They are further bounded at infinity, and coincide 
on the real axis z = f, where there is a single pole of order 
one located at z = 0. 

^Note that, under the assumption of finite second moment, 
K{f) 1 7-1 



lim 3? 



where the last equality straightforwardly follows by repeated appUcation of 
De L'Hopital rule, and from K' {0) = i2w and K"(0) = -47r2(l + 7). This 
ensures that W{f) is well-behaved. 
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An application of the analytic continuation theorem ||29l 
will allow to glue together the two functions in the upper and 
lower half planes, obtaining a function which is analytic in 
the whole plane, except for the single pole of order one at 
the origin. The generalized Liouville theorem f29 \ defines the 
only admissible form that such a function can assume; c/z, 
where c is a constant to be determinecjfl 

Restricting to the real-axis only, we finally get: 

^ [W{f)e-^-'ft - - 



f 



(18) 



Computing c is made possible by the condition that the 
sought u{t) should be a probability density function, which 
is equivalent to U{0) = 1, or V+{0) - V-{0) = 1 - 0(0). 
Using eqs. (fTsT l, we can write 



V+if) = [l-Kif)][W{f)e 



V~if) 

pi2irSf 



= [l-K{mW{f)e 



-iTrSf] 



—irrSfi 



f 

K{f) 



I 



that, evaluated at / = 0, yield T^+(0) - V""(0) = -lAirc, 



where we used lim/ 



i-K(f) 



i2n. The where 



condition ^+(0) - V""(0) = 1 - f7(0) will thus give 

1 - n{o) 

c = I . 

If we repeat the above development by multiplying eq. (flSl l 
by the complex exponential e^'^^^ , we get a similar result, 
namel>0 



(19) 



1 - ^(/) 

y-{f) , 



[VK(/)e 



i7ri5/] 



C 

7 



(20) 



Putting together eqs. (fTSl l and ( |20] |. along with the found 
value of the constant ( fT9l ), gives the system of equations 



l-^(/) 

l-K{f) 

v-{f) 

l-K{f) 



[W{f)t 



[w{f). 



[Win 



171-5/1 



Anf 

1 - n{o) 

inf 

1 - n{o) 

inf 

1 - n{o) 



(n) 
(Hi) 

(iv) 



1-K{f) ' ' inf 

Solving for V^{f)/[l — K{f)] in equations (i) and (ii) 



gives 



[Wif)e 



177(5/"! ^ e~^^^-^ 



[Wif)e- 



-i7r(5/"| + giTTiS/ 



Ssmc{5f) 



1 - fliO) 



^Actually, according to the generalized Liouville theorem, the overall 
function should be equal to cq + ci/z. On the other hand, we are looking 
for a solution U{f) in the class of the functions which vanish at infinity, 
implying cq = 0. 

*The structure of the equation is such that the same values of the constant 
c = i[l- f2(0)]/(47r) is obtained. 



(Using (Hi) and (iv) gives identical results.) 

The LHS of the above is equivalent to low-pass filtering of 
W{f), namely J W{i')S smc[6{f — u)]du, so that, by further 
using the explicit form ( fTSl l of W{f), and the properties of 
we get the desired claim. • 

For later use, note that at LHS and RHS of eq. (|6]l appear 
Fourier transforms of functions that vanish outside the range 
[—5/2, 6/2]. In view of the sampling theorem |30|, the samples 
taken at h/6, h integer, define the whole functions. These 
samples are 



n{h/s) 



1 - K{u) j 



(5sinc((5j/ — h)dv 

, 1 - m r 



where Ih — ^ for h — and Ih — otherwise. 

Also the unknown function ri(/) is bandlimited, so that 
it can be represented by the Shannon series = 
J2k ^{k / 6)smc{5 f ~ k). Substituting into the above equation 
we get 

^Ahk^{k/S) = -h, (21) 



k 



Aqo = 1 



K{v) 



4 



6s,\nc^ {5i')dv 



Akk = 1 



2 / SR { i5ff^) } 5^,mc^{6v - k)dv, k^O 

Ahk = 2 J 3? I i^K{i^) } '5sinc(5i^ — h)smc{Si' — k)dv, h ^ k 

(22) 

having used the orthogonality property of the sine functions. 

Thanks to the results of appendix |B] the above integrals in 
the Fourier domain can be expressed in the time domain as 
given in eqs. (l8]l-(fTTli. We are now in the position of proving 
the remaining claims. 

B. Proof of Theorem 3 

Let us consider only the term h = in ( |2TI ): 



AooniO) + Y,AoMk/6) = -. 



(23) 



The rationale behind the approximation of Theorem 3 amounts 
to neglect the cross-terms Aok, for fc 7^ 0, which can be easy 
understood by considering the two limiting regimes. 
Consider first (5 <C 1. By triangle inequality 



Ofc 



< 



m{t) dt^O 



where the approximation follows by m(0) — 0. 

As to (5 ^ 1, from a renewal theorem for interarrivals with 
finite second moment |24|, we know that 

lim [TO(i) -i] = (24) 

Thus, from (fTOl t we write, for k ^ 



i{t) - t - 



7-1 



cos(27rfct/(5) dt, 
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which follows by t CQs{2iTkt/ 5) — 0. Then, by triangle 
inequality 



Ofc 



2 

< - 

- 5 



m{t) - t - 



7-1 



dt « 0, 



where the last approximation is a consequence of the Cesaro 
mean theorem and eq. (l24l i. 

We note also that \Vl{k/5)\ < n{0) for all k ^ 0, and 
consistently we neglect all terms with k ^ in eq. ( |23]) . 
Solving for fl{0) is now possible: 



V. Examples 

The analytical expression of C provided by Theorem 3 
turns out to be quite accurate for virtually all the interarrival 
distributions used in our simulation studies, many of which are 
typical of network applications. Part of these extensive com- 
puter investigations are summarized in Sects. IV-AI and IV-BI 
In fact, to find an example where the refinements offered by 
Theorem 4 provide meaningful improvements over C, we need 
to choose carefully the kind of interarrival distributions, as 
discussed later 



2A, 



00 



m{t)dt 



yielding, in view of eq. (|5]l, the desired result. 



C. Proof of Corollary 1 

We know that the true embedding capacity tends to unity 
as 5 diverges. In order to quantify the convergence rate, we 
consider the Umiting behavior of 



1-C = 



l + f/>(i)dt-5 



1 + I /o Mt) dt 



(25) 



[t) dt^ - [m{t) - i] + (5 - 7 - 1 + (5, 



Now, 



5 J, 

by simple application of the Cesaro mean theorem and of the 
renewal theorem used before, see (l24l . From (|25| . we get the 
desired result: 

lim [1 - C]5 = 7. 

5— fC30 



D. Proof of Theorem 4 

Let us consider the series in eq. (l2ll . By truncating it to 
27V+1 terms, we get the following representation, with h,k 

[-N,N], 

N 



Ahk^(k/5)^-h 



Let A denote the matrix made of the Ahk'&- Recalling that we 
are not interested in computing the whole function r2(/), but 
just its value at the origin, we get 

It is possible to explicit the solution for = 1. Using the 
symmetries involved, we easily get 



{A-'}oo = 



1 



Ann + 2 ■ 



Aai-Aii 

which, along with eq. Q, yields the desired result (fTZt . 



A. Examples with Capacity in Closed Form 

We start with studying some well-known interarrival distri- 
butions, for which the renewal function is available in closed 
form. 

• Erlang family 

The Gamma random variable u-PDF is 



k{t) = e- 



r(0 



t > 0, ^ > 0. (26) 



When the parameter ^ is integer, the interarrival distribu- 
tion belongs to the Erlang family, a case for which the 
renewal function has been computed in closed form ||3T| 



«-i ah 



with e = . 

The (approximation of the) embedding capacity is ac- 
cordingly 

C 



?-l nh ( 



1 _ e-e(i-e'')AA^ 



Weibull distribution 

The u-PDF for the Weibull random variable is 

6-1 

t > 0, 6 > 



where cr = [r(l + 1/6)] ^. The pertinent u-RF is 

°° (-l)"-ia„ [r(l + l/6)i]"' 



(27) 



^(0 = E 



r(l + n6) 



where the coefficients a„ are defined recursively by 

n-l 

ai = ai . . . a„ = a„ - E '^3 o-n^j, 



with 

This yields 

C= 



r(i + ?i6) 



AA 



1 + 2E 



(-l)"-ia„ [r(l + l/&)AA]^ 



lb ■ 



(28) 



T{l + nb){nb+l) 
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Exponential 



Eriang ^=2 




XA 
Eriang ^=3 



XA 
Uniform 




Fig. 4. Examples of traffic models for which the renewal function admits 
simple closed forni. Dots refer to computer simulations of the embedding 
capacity and lines refer to analytical formulas. 




Fig. 5. Examples of different traffic models. Continuous curves refer to the 
approximation C in Theorem 3, eq. (|7), while dots are obtained by computer 
simulations. 



u-PDF is 



Uniform distribution 

The u-RF for uniform random variables can be obtained 
by iteratively solving the renewal equation 1.24.1 . yielding 



m{t) 



e*/2 -1 0<t<2 
et/2 _ 1 „ (I _ 1) e*/2-i 2 < t < 4 



with similar expressions for successive intervals of length 
2. The resulting capacity is 



C= < 



TS- 

4e~ -AA-4 



4e^+2e^"^(4-AA)-AA-8 



AA e [0, 2] 
AA e [2,4] 



In Fig. |4] the above expressions for the capacity are com- 
pared to numerical simulations. We see that the matching 
between the approximate analytical formulas and the results 
of computer simulations is excellent. 

B. Other Distributions 

Even when the renewal function is not known explicitly, 
there exist many numerical ways to compute that. Some 
methods exploit the definition of the renewal function in terms 
of interarrival distribution [24], other approaches are based on 
the interarrival density, and even others exploit the Fourier 
domain. 

For instance, let us consider again the Gamma family ( |26] l. 
and assume that ^ = 0.3. In this case, it is particularly 
convenient to use the Fourier domain expression for Aqo 
given in the first equation of ( |22] |. Computing numerically the 
involved integral, we get the capacity plotted in Fig. |5] 

A case of special interest for network applications due to 
its tail behavior is the Pareto interarrival distribution, whose 



m = 



b/{b-i) 



(29) 



Figure |5] shows the embedding capacity, again obtained by 
numerical integration of Aqo in < l22l l. for the Pareto distribution 
This distribution exhibits finite second moment whenever 
the shape parameter 6 > 2. Accordingly, in the numerical 
simulation, we first test the case 6 = 3, which falls in the class 
considered in the assumptions of our theorems, see Fig. |5] 
Then, we explore by simulation a case with infinite second 
moment, that is, 6 = 1.5, and Fig. |5] reveal that the accuracy 
of the formula is still excellent. 

In all the cases examined so far, there is no doubt that 
the expression C is quite accurate for any practical purposes. 
We would like to present an example in which the analytical 
formula (|7]i of Theorem 3 is less accurate and the following 
shifted exponential distribution offers this opportunity. Let us 
consider the following u-PDF for the interarrivals: kit) = 
e^T^, for t > a, and < a < 1. 

The embedding capacity is displayed, together with the 
simulated data, in Fig. |6l As it can be seen, the agreement is 
perfect in the range AA < a and is quite good for large AA; 
this is expected in view of Remark 2, and the arguments used 
in the proof of Theorem 3. However, for intermediate values 
of the product AA, the approximation C is not satisfying. 

Thus, we resort to the refined approximations offered by 
Theorem 4, and the results are shown again in Fig. |6] where 
the solutions obtained by using = 1 (that is, eq. (fT2]i). and 
N ~ 2 (this case being solved numerically), are displayed. 
As it can be seen, the partial inaccuracy of the approximation 
C is remediated with the adoption of eq. (fTSI i. Adding more 
terms (i.e., > 1) gives negligible improvements. 

VI. Ordering of Embedding Capacities 

In this section we show how the embedding capacity C 
can be used for comparing different renewal processes in 
terms of their embedding capabilities. Let Xi and X2 be two 
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Fig. 6. Example of a shifted exponential distribution, with a = 0.8. Dots 
are obtained by computer simulations, while continuous curves refer to the 
different analytical approximations for C* in Theorems 3 and 4. Specific ally , 
we display (i) C, (it) the linear system solution with A'^ = 1 of eq. )12K 
and [iii) the linear system solution for N = 2. The latter two curves are 
superimposed. 



non-negative random variables with the same average value 
E[Xi] = E[X2] = 1/A, and with cumulative distribution 
functions denoted by Fi{-) and F2{-), respectively. The fol- 
lowing definitions and results are classical in stochastic order 
literature, and can be found in ll33l . ||241 . 

Definition 7 (Variability or convex ordering) The random 

< „ X2, if 



variable Xi is less variable than X2, written Xi 
E[(j){Xi)] < E[(j){X2)] for all convex functions 



(30) 

o 



provided that the expectations exist. 

Known results L33J (Sufficient and Necessary Conditions 
for convex ordering) For non-negative random variables Xi 
and X2, with E[Xi] = £[^2] = 1/A, the condition Xi <y X2 
is equivalent to each of the following: 



Fi{t)dt > / F2{t)dt, for all x, (31) 
"'0 

Flip) > F2ip), for all [0,1]. (32) 

In the above, Fi(p) and F2{p) are the so-called Lorenz curves 
of the random variables Xi and X2 defined as 

rp 

Fi,2{p) = A 



' / Fi2{'U')du, for all p e [0, 1] 



Intuitively, Xi <„ X2 if Xi gives less weight to the 
extreme values with respect to X2- One way to get this is 
just to ensure that E[0(Xi)] < E[(/)(X2)] for convex 0, as 
stated in ( l30l l. That's why this kind of stochastic ordering 
is also known as convex ordering. It is also obvious that 
Xi <^ X2 ^ VAR[Xi] < YAR[X2], and hence Xi has a 
dispersion index smaller than or equal to that of X2, a fact 
that plays a major role in the regime of A 1/ A, as seen in 
Corollary 1. 

The following theorem formally relates the classical con- 
cept of variability ordering to the embedding capacity in a 



straightforward and intuitive way: less variable interarrivals 
yield a larger embedding capacity. 

Theorem 5 Let Ci and C2 be the approximate embedding 
capacities for i.i.d. renewal processes with interarrival distri- 
bution Xi and X2, respectively. Then 



Xi<,X2 ^Ci>C2. 



Proof. The u-RF's of Xi and X2 can be represented as 







mi 



00 00 
it) = {^^f ^ < 4 ' W = E Pr { A5f ) < t} 

i=l 1=1 

(33) 

Let us focus on the single terms of the series. By assumption 
Xi<yX2, implying, in view of (|3T| i. 



AA 



Pr 



{A5«<i} 



dt < 



AA 



Pr 



dt. 



Since the variability ordering is closed under convolution (see 

e.g., m) 



AA 



Pr 



{xsip<t} 



■ AA 

dt < I Pr 





{xsi^^<t} 



dt. 



This implies 

r±Pr{xsl'^<t} dt< r±Pr{xsl'^<t} dt. 
Jo Jo 

Applying Beppo Levi's monotone convergence theorem |34), 
we are legitimate to exchange integration and limit, yielding 

AA i-XA 

mi{t) dt < / 1712 (t) dt 
Jo 

which, in the light of eq. (|33] |. gives Ci > C2. • 

A. Ordering w.r.t. Poisson 

It is of special interest to compare the given renewal process 
to Poisson traffic, and this can be conveniently done by means 
of our analytical approximation. To do so, let us define two 
special categories of interarrival distributions. 

Definition 8 (NBUE/NWUE classes) A non-negative random 
variable X is called New Better than Used in Expectation 
(NBUE) or New Worse than Used in Expectation (NWUE) 
ifm: 

NBUE E[X - s\X > s] < E[X] V.s > 0, 
NWUE E[X - s\X > s]> E[X] V.S > 0. 

Due to the absence of memory, the exponential distribution 
is such that E[X - s|X > s] E[X], and it belongs to both 
classes. 

Corollary 2 (Capacity Ordering in NBUE/NWUE classes) 
Let C[.^gij£, Cf^-\j\/u£, and Cexp denote the embedding capac- 
ities given by (0 for interarrivals from the NBUE class, the 
NWUE class, and the exponential distribution. The following 
relationship holds: 

^NWUE < '^exp < Cf^BUE- 
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Proof. Thanks to Proposition 9.6.1 in [ED, the NBUE (resp. 
NWUE) distributions can be shown to be less (resp. more) 
variable than the exponential, implying the claimed result as 
a direct consequence of Theorem 5. • 





u-PDF k{t) 


Ordering Relationship 


GAMMA 


Eq. (26) 


?i > 6 ^ C-i > Ca. 


WEIBULL 


Eq. (27) 


fel > 62 ^ Ci > C2. 


PARETO 


Eq. (29) 


6l > fe2 ^ Ci > C2. 


LOGNORMAL 


Eq. (34) 


(71 < 0-2 ^ Cl > C2- 



TABLE I 

Summary of the relationships between classical and 
embedding capacity ordering for typical distributions. 



O 
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CO 
O 



■a 
■D 

_g 0.3 

E 

UJ 

0.2 
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□ Scrambled data 
Weibull (ttieor.) 



O o 



o <> 



absolute errors (different tranches) ^ 



?i A 



B. Ordering within the same distribution class 

The relationship between embedding capacity ordering and 
classical ordering of random variables allows easy comparison 
of distributions within the same class. 

For Gamma and Weibull random variables, it has been 
shown that the Lorenz curves are monotonically increasing 
with the shape parameters ll35l . Thus, larger shape parameters 
give higher embedding capacities. It is also easy to evaluate 
the Lorenz curve of the Pareto random variable with a u-PDF 
given in ( |29] l 

L{p) = h{l-p) [1-{\-pY-^'']+p, 

as well as that of the Lognormal random variable 

L(p) = $($-i(p)_a), 

whose u-PDF is 
1 



k{t) = 



cxp 



2a2 



t > 0. 
(34) 



Both functions exhibit monotonic behavior with respect to b 
and CT, respectively. 

Therefore, using eq. (l32t allows easy (convex) ordering of 
the interarrivals, which in turns induces an ordering of the 
embedding capacities thanks to Theorem 5. The results are 
summarized in Table Jl 

VII. Experiments with Real Network Traces 

In this section, we present some numerical tests run on 
real traces. Specifically, we downloaded the TCP packet ar- 
riving times (traces lbl-tcp-3.tcp and lbl-pkt-4.tcp) gathered 
at the Lawrence Berkeley Laboratory, that were originally 
used in [36]. Following |36| and [4], we extract packets 
corresponding to Telnet connections (obtained from hearing 
communication on port 23). The pipeline for the real data 
processing is as follows. 

• The inspected traces correspond to traffic patterns col- 
lected in two different days. We consider the aggregate 
traffic, that is, we do not extract informations pertaining 
to the single hosts. 

• We emulate the scenario of two mutually independent 
point processes, by using two tranches of 10"* packets 



Fig. 7. Embedding capacity curee of Telnet data, for a pair of tranches 
selected as described in the main text. In the lower part of the plot, the 
absolute error between empirical and theoretical capacity is displayed, for a 
broader set of different tranches. 



each, extracted from traces lbl-tcp-3.tcp (source node) 
and lbl-pkt-4.tcp (relay node). 

• By means of a moving average filter over 10^ packets, 
we select over the two traces candidate tranches having 
comparable rate^. Without loss of generality, we scale 
the data by dividing the interarrivals by the sample mean 
computed over the union of the two tranches. 

• We run the BGM algorithm on the selected tranches, with 
fixed (dimensionless) observation time t — 9000. 

« We also run the BGM algorithm after scrambling the 
interarrivals, in order to remove statistical dependencies 
between them, namely, to enforce the renewal assump- 
tion. This is purely for testing the accuracy of the found 
formulas. 

In order to compute theoretical capacities, we need a candi- 
date distribution for the interarrivals. We accordingly fit the 
empirical interarrival CDF of each tranche, and find that the 
Weibull distribution works generally well, that is perhaps not 
unexpected, see, e.g., ll37l and ll38l . 

Consider now the capacity curves in Fig. |7] The exper- 
imental curves for capacity refer to one pair of tranches 
where the Weibull fit is accurate, and the two empirical CDFs 
are close to each other, complying with the assumption of 
identical distribution across nodes. The theoretical curve is 
drawn by (l28l l, where the shape parameter b is computed over 
the union of the two tranches. 

For the scrambled data the theoretical approximation C is 
excellent. As to the (non-scrambled) real data, a first evidence 
is that, up to values of AA in the order of unit, the curve 
matches the theoretical approximation well. On the other hand, 
a discrepancy emerges at larger values of the product AA, due 
to possible dependencies among the interarrivals. 

^With this selection procedure, the tranches extracted from a given trace 
might also overlap. Obviously, this does not alter our analysis, in that we 
only need independence between the source tranche (extracted from trace 
lbl-tcp-3.tcp), and the relay tranche (extracted from the independent trace Ibl- 
pkt-4.tcp). 
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A more complete picture is obtained by applying the above 
procedure to different tranches, irrespectively of the goodness 
of the Weibull fit, and of the similarity between the empirical 
distributions at the two nodes. The results of this latter 
analysis are summarized in the bottom part of Fig. |7] where 
the absolute error between the theoretical formula and the 
empirical capacity is displayed, only for the case of real data. 
(Again, scrambling reduces the error, this is not shown in 
the plot.) The points marked with darker circles refer to the 
tranche pair used for drawing the capacities displayed in the 
uppermost part of the plot. As it can be seen, the theoretical 
approximation follows the empirical capacity closely at small 
AA. Also in this case, a discrepancy is observed for moderately 
large values of AA, with an absolute error staying in the order 
of 10-1. 

Summarizing, a main behavior seems to emerge — that 
the theoretical predictions are very accurate for real data 
well modeled by renewal processes, corroborating the whole 
theoretical machinery for embedding capacity computation, 
and that the possible statistical dependence among packet 
interanivals can be neglected for tight delay constraints, up 
to delay values in the order of the mean interarrival time. 

VIII. Conclusion 

We consider the problem of matching two independent 
and identically distributed renewal processes, according to a 
bounded delay criterion, with applications to communication 
network scenarios. We introduce the concept of embedding 
capacity, and provide fully analytical tools and approximations 
to evaluate it, relying upon the Riemann-Hilbert theory. An 
exact evaluation of the capacity is reduced to a manageable 
integral equation, that can be solved to any degree of ap- 
proximation by inverting a highly structured linear system. 
The main finding, however, is a simple approximated formula 
of the embedding capacity that involves the renewal function 
of the underlying processes. The approximation is excellent 
for virtually all the cases of practical interest that we have 
investigated, part of which are reported in the paper Even 
when this is not strictly true, we provide closed-form solutions 
for first-order correction. 

The analytical formula highlights the role played by dif- 
ferent renewal parameters: for large AA only the dispersion 
index matters, while embedding capacity ordering is induced 
by the stochastic variability of the underlying interarrivals. 

The experimental analysis carried on real network traces 
reveals that the accuracy of the analytical expression is good 
for tight delay constraints, up to AA in the order of unit. For 
larger delays, a partial inaccuracy is seen, and we show that 
this should be ascribed to statistical dependencies unavoidably 
present in real traffic patterns: the renewal model is failing, 
rather than the proposed analytical approximation. 

The abstract concept of matching between point processes 
arises in a very large number of contexts, and we feel that our 
findings can represent a contribution to these fields. To broaden 
further the horizon of potential applications, refinements and 
improvements of the approach can be considered. These the 
case of different renewal processes at the two nodes, the 



extension to non-renewal point processes, to multi-hop flows, 
and to the case of multiple input/multiple output relays, see lO, 

Appendix A 
Proof of Theorem 1 

We first justify the embedding capacity formula dU. Assume 
for now that the frequency for Z„ to fall inside the interval 
[0, A] converges a.s. to a constant p. Then since each Z„ 
outside [0, A] represents a chaff point whereas each Z„ inside 
the interval represents a pair of flow points, we see that the 
fraction of flow points embedded by BGM converges a.s., and 
the limit, i.e., the embedding capacity, is given by 2p/{l +p). 

On the other hand, by Theorem 17.1.7 in ||39], if {Z„} 
is a positive Harris recurrent Markov chain, then i) p ^ 
lim„^oo ^I]"=o-^[o, A](^j ) exists a.s., where /[o, a](2^) is 
the indicator function, and ii) p can be computed from the 
invariant PDF h{t) hy p = h{t)dt. By definition, if h{t) is 
the solution to eq. (|3]l, it will be invariant under the transition 
dU, i.e., it is an invariant measure. The positive Harris property 
of the chain implies that h{t) is unique and finite, and thus 
can be normalized into a probability measure. It remains to 
prove the property of positive Harris recurrence. 

First, we show that the Markov chain {Z„} is V'- 
irreducible [39 1 (all the sets mentioned in the sequel are Borel). 
The assumption that BGM can match one pair almost surely 
implies that the interval [0, A] is accessible from any state 
almost surely, say L{z, [0, A]) = 1 Vz [39J. This rules out the 
cases where the asymptotic fraction of matched points depends 
on the initial state (where embedding capacity does not exist) 
and those where the embedding capacity is trivially zero. 

Let If be the Lebesgue measure constrained to [0, A], i.e. 
(p{A) = fi{An [0, A]), where ^ is the Lebesgue measure. 
Given PDF f{t), there must exist eo > such that f{t) > 6o 
for all t within some interval [to, to + eo], and thus 

m = / /(r)/(r - t)dT > S^eo - \t\) > 5l{eo - ei) 
Jo 

for all t e [—61, ei], where ei is a constant in (0, eg). Let 
5i := (5Q(eo— ei). Partition [0, A] into m [2A/ei] segments 
of length ei/2, as illustrated in Fig. [8] such that the transition 
density from any z S [0, A] to any point in an adjacent 
segment is greater than Si. For any set C with (p{C) > 0, 
let £2 be the Lebesgue measure of the minimum intersection 
between C and the ^-segments. Let z be an arbitrary point 
in [0, A] that is n segments away from C (n < m — 1) and 
Ii (i — 1, . . . ,n) be the ith segment from z to C, where I„ 
intersects with C. The n-step transition satisfies 

P'^{z,C)> / !o{xi-z)dxi / fo{x2- xi)dx2 - ■ ■ 
Jxi J12 

/ fo{Xn - Xn-l)dXn 

Jx^nc 

> (<5i^J 5162 >0. (35) 

This implies L{z,C) > for all z e [0, A]. Moreover, since 
L{z, [0, A]) = 1 for aU z, we have L{z,C) > for aU z. That 
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Fig. 8. Access C from z by hoping tlirougli ^-segments in [0, A]. 



is, any set with positive measure is accessible from anywhere 
within the state space with positive probability, implying 
that the chain is (/^-irreducible and hence T/i-irreducible for a 
maximal irreducibility measure according to [39J. 

Second, we show that {Z„} is Harris recurrent. Since it 
is -i/j-irreducible and [0, A]) > for all z, by Theorem 
5.2.2 in |39|, there exist > 1, a nontrivial measure V}^, and 
a nontrivial set C\ C [0, A] such that C\ is z^^-small, and 
hence z^^^-petite. For sampling distribution a{i) = 1/m (i = 
1, . . . , m), the transition kernel of the sampled chain from any 
z € [0, A] satisfies 

Ka{z,Ci) > -P"(z,Ci) > - (6i'-^)7il2, (36) 
TO m \ 2 / 

where we apply ( l35T l for C = Ci. Since n < m — 1, 
Ka{z,Ci) > i ((5iei/2)™ Sie2, independent of z for z E 
[0,A]. Therefore, Ci is uniformly accessible using a from 
[0,A]. By Proposition 5.5.4 in we prove that [0, A] 

is ;^a*(5fc -petite. The fact that a petite set [0, A] satisfies 
L{z, [0, A]) = 1 for all z for a ■(/'-irreducible chain implies 
Harris recurrence in the light of Proposition 9.1.7 in f|39l . 

Finally, we show positivity by drift analysis. Define the 
function 



V{z) = 2A ■ 



z - 
0, 



A, 



if z > A, 
if < z < A, 
if z < 0, 



where 1/ A is the mean interarrival time, and consider the mean 
drift defined in 1391 as 



dV{z) 



Piz,dy)V{y)-V{z), 



where P{z,dy) is the transition kernel of the chain. Define 
a set C2 = [—zq, a + zq] for zq sufficiently large such that 
fit)tdt - J^^^ f{t)tdt > 1/(2A). For any z > A + zo 
we have, after some sti-aightforward manipulations. 



dV{z) 



< 



-2A 


- fim- 


/ f{t)tdt + (z - 
Jo 


-2A 


(\f{t)tdt~ 

Jq 



f{t)dt 



z-A 
f(t)tdt 



zo+A 



< -1. 



The same holds for z < — zq. It is easy to see that, inside the 
set C2, dV{z) can be bounded by a constant, such that we can 
write 



dV{z) < -l + 6/c,(z), 



(37) 



with a suitable choice of h. Since the petite set [0, A] is 
uniformly accessible from C2, we can conclude that C2 is 
petite, and eq. (|37] i coincides with the drift condition {iv) 
of Theorem 13.0.1 in |f39l . whence, further observing that 
aperiodicity holds, we conclude that {Z„} is positive Harris. • 

Appendix B 
Linear system coefficients 

Let us introduce the so-called renewal density associated 
to the renewal function m{t), that is p{t) = dm{t)/dt. It is 
convenient to consider a symmetric version thereof, namely 
p{t) — p{t) + p{—t). It holds true that the Fourier transform 
of p{t) - 1 is given by 2 5r|j^J^}, see HO), BTl . 

Let us first consider the term Aqq in eq. ( |22] |. In view of 
Parseval's formula |42|: 

gfj I ^''"l I ^sinc^ (Siy)diy 
1 - K{i^) J 

[p{t)-l]{l-\t\/5)dt = 2 / p{t){l-t/S)dt-S, 
i-s Jo 

where we simply notice that the Fourier transform of the 

triangular window of width 26 is Ssinc^ (Sf). Integration by 

parts then gives 2 pit){l — t/S)dt = | m{t)dt, or 



^00 = 1 - + T 



m(t)dt. 



2 ■ S 

As to the evaluation of A^k in eq. (l22T i. fc 7^ 0, it suffices 
to use the shift property of the Fourier transform, yielding 

KM 



1 - K{iy) 



(5sinc (6iy — k)dv 



[p{t) - 1](1 - \t\/5) cos{27Tkt/S)dt 



s 



= 2 / p{t){l -t/ 6) 008(2^/ 6)dt, 
Jo 

that integrated by parts gives 
2 

A-kk = 1 + T / m{t)[cos{2TTkt/6)dt 
Jo 

+ 27rfc (1 - t/S) sin(27rfct/(5)] dt. 

Finally focusing on the terms Ahk in eq. ( |22] |. /i 7^ fc, it 
suffices to consider the even part of (5sinc((5/ — /i)sinc((5/ — fc), 
whose inverse Fourier transform is 



1 



f.5/2 



I ^-^2rr{h-k)f^- 
I J-S/2 
1/2 



f ^ I dT 



cos[27r(/i - k)T + 2Tikt/5] n(r - t/5) dr, 

-1/2 

n(t) being the rectangular window of width 1. The integral is 
zero for \t\ > 6. For t e (0,(5) we have 

nl/2 

cos[27r(/i - k)T + 2Trkt/S]dT 

it/s-1/2 

(^_l-^h-k 



2TT{h-k) 



' sm{2Trkt/S) - sin{2Tr ht / S) ] 



°This can be easily shown with the same technique used to prove uniform 
accessibility of Ci from [0, A]. 
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This gives 

Ahk = ,_^ 2 / p{t)[sm{2TTkt/S) - sm{2TTht/6) ] dt 



2'!r{h - k) 7o 

m{t)[h cos{2TTht/S) — k cos{2TTkt/6)]dt, 



{-l)h'k 2 rS 



(h-k) SJo 

where the latter is obtained integrating by parts. Equation ( fTOl i 
now follows as a special case, whence eq. (fTTI) is true. 
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